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Standard variational methods tend to obtain upper bounds on the ground state energy of quantum 
many-body systems. Here we study a complementary method that determines lower bounds on 
the ground state energy in a systematic fashion, scales polynomially in the system size and gives 
direct access to correlation functions. This is achieved by relaxing the positivity constraint on the 
density matrix and replacing it by positivity constraints on moment matrices, thus yielding a semi- 
definite programme. Further, the number of free parameters in the optimization problem can be 
reduced dramatically under the assumption of translational invariance. A novel numerical approach, 
principally a combination of a projected gradient algorithm with Dykstra's algorithm, for solving 
the optimization problem in a memory-efficient manner is presented and a proof of convergence for 
this iterative method is given. Numerical experiments that determine lower bounds on the ground 
state energies for the Ising and Heisenberg Hamiltonians confirm that the approach can be applied 
to large systems, especially under the assumption of translational invariance. 
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I. BACKGROUND AND MOTIVATION 

The determination of the ground state properties of strongly interacting quantum systems is of considerable interest. 
As only a very few models are exactly solvable, numerical approximation methods are of significant importance. A 
key challenge in this respect is the 'curse of dimensionality', namely the exponential growth of the Hilbert space 
dimension with the number of sites. Ingenious methods such as the density matrix renormalization group (DMRG) 
approach [2] can overcome this problem as they amount to formulating a variational problem over a well-chosen 
class of quantum states that can be parametrized efficiently i.e. with a polynomial number of parameters [3]. As the 
variation is over a subclass of all possible quantum states, variational methods of this type provide upper bounds on 
the ground state energy. Assessing the quality of the so obtained upper bounds on the ground state energy is then 
difficult. Hence there is considerable interest in the development of efficient methods for the determination of lower 
bounds on the ground state energy. One straightforward approach in this direction consists in the determination of 
Anderson bounds 0], i.e. the decomposition of a Hamiltonian H = J^ fc Hk and the subsequent determination of the 
ground state energy of each Hk whose sum then yields a lower bound on the ground state energy of H. While useful, 
this approach does not scale well with the system size as it remains exponential in the size of the support of the Hk- 
The application of DMRG to the Hk may overcome the scaling issue but amounts to finding an 'upper bound to a 
lower bound' and does not resolve the principal issue with DMRG. 

As early as 1940, however, it was suggested to avoid the use of wave functions altogether and rather consider 
reduced density matrices together with conditions that these arise from a valid quantum state [5] - the so-called 
N-representability problem [5J [7j . Varying over the reduced density matrices that are compatible with a valid physical 
state then allows for the determination of the ground state energy. The exact ground state energy can be obtained 
if all conditions for the N-representability problem are known. If only partial conditions are specified, then one 
minimizes over too large a set of reduced density matrices and one obtains a lower bound on the ground state energy. 
The analytical treatment of this problem is highly challenging and at the time numerical solutions were hard to come 
by. Further progress was made by Mazziotti [5J 15] , who pointed out that the problem can be formulated as a semi- 
definite programme and use can be made of the guaranteed and certified convergence of the primal-dual interior-point 
algorithm for semi-definite programming [101 lll| . A general mathematical framework for this type of optimization 
problems and applications was given in |12j . Unfortunately, the primal-dual interior-point algorithm and related 
methods do come with a considerable overhead in computation time and, crucially, memory that is making it hard 
to scale to very large systems. Lately, new algorithms to address the optimization problems of finding lower bounds 
on ground state energies have been developed [T51 H"4"] and considerable improvements towards the interior-point 
algorithms have been demonstrated. 

In this work, we are exploring further the idea of relaxations of the ground state energy finding problem that 
allow us to determine lower bounds on the ground state energy in an efficient manner. To this end we formulate 
the problem, introduce a number of simplifications including the explicit treatment of translation invariance and 
symmetries and point out additional potential routes for numerical simplifications. In addition, we present a novel 
algorithm, principally a combination of a projected gradient algorithm with Dykstra's algorithm |15H17j , for solving the 
resulting optimization problem in a memory-efficient manner and prove its convergence. Then, in order to demonstrate 
the capabilities of this approach, we present numerical experiments that determine lower bounds on the ground state 
energies for two one-dimensional (ID) lattice models, namely the Ising and Heisenberg Hamiltonian, which confirm 
that the approach can be applied for large systems especially when implementing translational invariance. Note 
that although we restrict our numerical analysis to ID lattice models it has been demonstrated that the variational 
determination of ground state energies considering only constraints on reduced density matrices applies to higher 
spatial dimensions [18] and to non-periodic systems such as molecules in quantum chemistry |13) as well. We complete 
this work with a discussion and outlook. 

II. BASICS 

We begin by explaining the general ideas underlying the method for the determination of lower bounds on ground 
state energies following [5] and then specialize to the case of fermionic systems, for which later sections will provide 
numerical examples. Analogous ideas can also be implemented for bosonic systems or general spin-systems. 

A. Generalities 

The determination of the ground state energy Eq of a many-body system composed of N particles and described 
by the Hamilton operator H can be formulated as the minimization problem 



E = min {tr [Hp] \ p > 0, tr [p] = 1} 



(1) 



3 



This approach is pushed to its limits rather rapidly as it uses the positivity constraint on p, which grows dramatically 
with the system size, e.g. for particles with spin-1/2 or spinless fermions p is a 2 N x 2^ matrix. Hence, this problem 
cannot be solved efficiently due to the exponentially increasing parameter space. One technique to deal with this 
difficulty is to parametrize a specific family of candidates \<j> (d?)) for the ground state by a set of parameters a. One 
then minimizes the expectation value of the Hamiltonian with respect to this family of states. As this family of states 
does not encompass the entire state space, this method will find a state \<f>o) which is optimal in {\<f> (a))} s and which 
yields an upper bound on the true ground state energy. This is the route taken by a variety of variational methods 
including, for example, DMRG [TH3]. 

To overcome the scaling issues and to provide lower rather than upper bounds on the ground state energy, we need 
to relax the positivity constraint p > 0, replacing it by a weaker set of constraints which scale only polynomially in 
the system's size [51 [T3] . This enlarges the class of states over which one varies the energy and hence yields lower 
bounds on the ground state energy. To this end we use the fact that 

p > tr [O f Op] > (2) 

for all operators O. Demanding only 

(O^O) = tr [0*Op] > (3) 

for some operator O does not, in general, imply p > 0. For some subset of operators {Ok}k, we can then define 
O = J2k a kOk and the requirement that (O'O) > for all choices of a k is then equivalent to the positivity of the 
matrix X with components X^i — (OtO;), imposing a particular N-represent ability condition. Hence, replacing in 
the original optimization problem eq. ([!]) the condition p > with the weaker condition X > and omitting the 
normalization tr[p] — 1, we find that 

E > min {tr [Hp] \ X > 0} . (4) 

The present formulation still contains the density operator p explicitly in the computation of tr [Hp]. Judicious choice 
of the set {Ok}k alleviates this shortcoming. Indeed, writing H = ^^hkiHi-i, where hki € C are the coefficients 
of the operators Hki in the Hamiltonian H, one ensures that the set {O k } k contains all operators such that for all 
k,l the expectation value (Hki) 1S m {{0\Ol)}kl- Note, furthermore, that the entries of X will generally not be 
independent but suffer some linear equality constraints, for example due to (anti)-commutation relations between 
operators. Denoting by C the set of all Hermitian matrices that satisfy X > and any other linear constraints, the 
energy minimization problem can be written as [8, 9 J 

E ° > min j^/ifciXfe, | Xecj, (5) 

which is evidently a semi-definite programme [TT and efficiently solvable as long as the set {Ok}k contains a number 
of operators that scale polynomially in system size. The decomposition of H into Hki , and hence the choice of the set 
{Ofcjfe, is not unique and the optimal choice is not obvious. From a practical point of view however, the systematic 
choice presented in the remainder of the paper leading to the so-called p-positivity conditions [8, 9] does perhaps 
provide the most straightforward formulation. 



B. Bosons and Fermions 



For concreteness, let us now consider the class of bosonic or fermionic operators defined by 

O = ^ S ><ak + ^ ffcaj, + a k ia k ai + ^ Pkia\ai + ^ lkia k a v ( 6 ) 

where 6k, £fc, cau, Pu, "f k i € C are arbitrary complex variables. The N-representability conditions arising from this class 
of operators are called 1-positivity and 2-positivity conditions, which is a terminology used to reveal the powers of the 
creation and annihilation operators a\ and in O. In general, allowing polynomials of order p in the second-quantized 
operators leads to restrictions called p-positivity conditions [H [S] ■ For the particular choice of operators of the order 
of 2, requiring the positivity of (O^O) is equivalent to 
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where the N x N matrices comprising the second moments are defined as follows: 

Tk,i = (Ofcffli)i 
Sk.i = (fflfcof), 

U k .i = (a k ai) (8) 
and for the third and fourth moments, we find the N 2 x N and N 2 x N 2 matrices 



(9) 
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In addition to the positivity constraint on the 27V + 3 TV 2 x 2N + 3N 2 Hermitian matrix X E %, where H denotes the 
vector space of all Hermitian matrices with appropriate dimension, there are linear constraints relating its entries. In 
fact, the fermionic anti-commutator relations 

{a k ,ai} = |<4, a|| = 0, ja fc , a[ j = <5 M (10) 

or the bosonic commutator relations 



[a k ,ai]= a{,aj =0, a k ,a] =4,; (11) 

force some entries of X to be equal, to differ by sign or to be a linear function of others. 

Now, let C be the set of all positive semi-definite Hermitian matrices of the form eq. |7| and satisfying the 
additional requirement that the entries of X obey the constraints given by the commutation or anti-commutation 
relations. Further, let E (X) : W — > K be a function expressing the expectation value of a Hamiltonian with up to 
4-mode terms, i.e. E (X) = (H). Then, the original optimization problem eq. ([!]) can be replaced by 

Eo > min{£pO \X E C] , (12) 

where the relaxed optimization problem scales polynomially in the system size TV. Note that it is straightforward to 
extend this formulation to higher moments, e.g. 3-positivity conditions. In the following, we are going to present 
scenarios and strategies where the number of free variables can be decreased considerably. 



C. Consequences of Superselection Rules and Particle Number Conservation for Fermions 



For the remainder we will consider fermionic systems. In that case, superselection rules require that the Hamiltonian 
commutes with the parity operator given by P — ~ 2aJ.a/ c ). As a consequence, only Hamiltonians comprising 

terms with an even number of creation and annihilation operators are physically allowed. This in turn ensures that 
only expectation values of products of an even number of creation and annihilation operators are non-vanishing. 
Hence, we find that we can restrict X to take the form 
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(13) 



and to satisfy the constraints dictated by the fermionic anti-commutator relations. See appendix [B] for a list of 
constraints on the entries for a formulation for up to fourth moments in the fermionic case. 

If the Hamiltonian conserves the particle number, a further simplification ensues. In this case, each term in the 
Hamiltonian will consist of an equal number of annihilation and creation operators. Hence, the Hamiltonian is 
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they concern operators with the same number of annihilation and creation operators so that X can be assumed to 
take the form 



X = 
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(14) 



Note that the elements of the matrices T,S are coefficients of the one-electron reduced density matrix (1-RDM), 
while the matrices M, R, Q are representations of the two-electron reduced density matrix (2-RDM) [5] . Restricting 
matrices of the form eq. ( 14 ) to be positive semi-definite is equivalent with the positivity of the individual matrices 



T, S and M, R, Q. This constraint on the matrices containing the second moments defines the 1-positivity conditions, 
while the positivity of the matrices containing the fourth moments generates the 2-positivity conditions 8, 9 . For 
concreteness we list the matrices comprising the constraints on a system composed of two spinless fermions and the 
additional assumption of particle number conservation in appendix [El 



D. Exploiting Translational Invariance 



So far we considered optimization problems for up to 2-positivity conditions which reduced the number of free 
variables to scale as N 4 . Apart from superselection rules and the particle number conservation, no specific assumption 
on the structure of the underlying Hamiltonian has been exploited such that in principle this procedure is applicable 
to a wide range of different systems, e.g. quantum chemistry calculations [13] or higher spacial dimensions [T5]. Since 
our intention is to apply this method to large lattice models, we now discuss a further simplification of the system 
to decrease the number of parameters, namely the translational invariance. Translation invariance with periodic 
boundary conditions of the physical system reduces the number of free variables considerably. In fact for the second 
moments, i.e. for the matrices T, S and U, this can be exploited immediately as these matrices then obey Ti j = U—j, 
etc. Consequently, for translational invariant systems with periodic boundary conditions, the matrices comprising the 
second moments become circulant matrices [19] and can therefore be diagonalized by a discrete Fourier transform 

Vkl = J_ e -2ni(k-l)(l-l)/N (15) 



for k, I = 1, . . . , N. Reducing the number of variables for the fourth moments is not so transparent as for the second 
moments, but manageable. To do so, note that the ordering of the entries of the matrices M, G, H, R, I and Q is not 
fixed but can be chosen arbitrarily. Remember that the indices of these matrices are of the form [(k, I) , (to, n)] for 
k,l,m,n = 1,...,N. Rearranging the two-digit indices labeling the rows and columns in the following way: 

(k,l) ={(1,1), (2,2),..., (N,N), 

(1.2) , (2,3),..., (TV, 1), 

(1.3) ,(2,4),...,(7V,2), (16) 

(1,N),(2,1),...,(N,N-1)}, 

the matrices comprising the fourth moments decompose in a block structure where each block can be diagonalized by 
the unitary matrix V^j given above. For example, we find for the matrix M 

M = UD M U\ (17) 



where U = diag (V, V, . . . , V) with V given by eq. ( 15 ) and 



D M = 
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^£)N,l ]jN,2 ]jN,N 



G 



with diagonal matrices D l ' J , i, j = 1, . . . , N. Hence, for translational invariant systems the matrices contained in set 
C will have the following form: 
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where the matrices denoting the second moments in the upper left block are diagonal matrices and the matrices 
considering the fourth moments have a block diagonal structure as described above. These considerations reduce the 
complexity of the constraints considerably and hence make the numerical study feasible for large systems. 



E. Varying System Size for Higher Moments 



As the treatment of the fourth-order correlations is computationally costly, it may in some cases be advantageous 
to consider fourth-order constraints on a smaller lattice, while the computationally cheaper second-order constraints 
are treated on the full lattice. In particular, while T, U and S are N x N matrices, M, G, H, R, I and Q are N 2 x iV 2 
matrices leading to the unfavourable scaling. To overcome this issue, it might therefore be beneficial to restrict higher 
moments to smaller subsystems; that is, for the second moments one would consider system size iVjj, while for fourth 
moments one would consider N4 where N2 N4. Reducing the support of the matrices M,G, H, R, I and Q, on 
the one hand, will lead to disproportionate savings in memory and computation time and, on the other hand, will 
certainly lead to a decrease of the lower bounds. 



III. THE PROJECTED GRADIENT ALGORITHM 

So far we have described how to relax the original optimization problem eq. ([I]) and discussed how to include 
symmetries to reduce the number of free variables. Now we are going to describe an iterative algorithm for addressing 
this minimization problem. Recall the minimization problem we have to face, that is, 

minimize E (X) 

subject to X eC ' 

where again the function E (X) : H — > K expresses (H) and the set C consists of all positive semi-definite Hermitian 
matrices X fulfilling the constraints dictated by the (anti)-commutator relations and potential other symmetries of 
the system. Note that the function E (X) is indeed linear in the elements of matrix X and hence could be written as 

E {X) = tr [G ■ X] + c, (21) 

where G is a Hermitian matrix and c G K. Furthermore, the set C is the intersection of the convex cone consisting 
of positive semi-definite matrices and the affine set defined by the linear constraints on the entries. Hence C is 
convex, too [11] . It is well known that this optimization problem is in fact a semi-definite programme and as such 
can be solved efficiently. Standard primal-dual interior point methods for semi-definite programming come with a 
considerable overhead in memory and time which makes their application to large systems challenging. Recently, 
two algorithms addressing the problem of ground state energy estimation were developed which scale better than the 
conventional methods [T31[T3]. Here, we continue the research for efficient algorithms and formulate an alternative 
numerical method whose efficacy we demonstrate in later sections. 

To address the issues described above, we suggest to solve the minimization problem by a projected gradient 
algorithm [2D]. Since this algorithm relies heavily on projections onto convex sets, we need to define the latter. Let 
H be a vector space and A E H be a convex set. The projection of X g H onto A is defined by 

P A (X) = a.rgmin{\\X ~Y\\ F \Y e A}, (22) 

where || • ||f : W — > R denotes the Frobenius or Hilbert-Schmidt norm. As the name suggests, the projection determines 
the nearest point in the convex set A with respect to a given norm. With this, we claim that the solution of problem 



eq. ( 20 ) could be found by iterating 



Xk+i = Pc {Xk - olG) , 



(23) 



E* 



FIG. 1. Illustration of the gradient projection algorithm 1 23 1. The convex sets A and B have a nonempty intersection, which 
is shown shaded. The straight parallel lines clarify the level sets of the function E (X) such that E* < E' < E" , where E* is 
the minimum value the function attains on the intersection. The projections are represented by arrows to demonstrate their 
directions. Note that these projections are computed using Dykstra's algorithm, which itself uses the projections onto A and 
B. The algorithm is initialized with the matrix Xq. 



where the algorithm is initialized by a Hermitian matrix X G % and a > is an arbitrary real and positive number. 
Furthermore, Pc ■ T~i — >• C denotes the projection onto set C, i.e. the intersection of the positive semi-definite Hermitian 
matrices and the set given by the linear constraints. The scheme of this algorithm is visualized in figure [I] 
We find that the following lemma holds. 

Lemma 1 Let E : H — >• R be a linear function of the form E (X) = (G, X) + c where G G %™ xn and c G M is a 
constant. Assume that the function E (X) is lower bounded on the closed convex set C, i.e. there is an X G C such 
that E(X) < E (X) for all X G C. Let E k = E (X k ) be the sequence determined by the iterations of the algorithm 



eq. (231. Then E k =E(X k ) -> E* where E* < E (X) for all X G C. 



A detailed proof of this lemma is presented in appendix [A] 

So far we have shifted the problem of minimizing the function E (X) with respect to the convex set C to the 
problem of computing the projection onto C. As pointed out above, this set is the intersection of two convex sets. 
This complicates the computation of the projection because in principle one has to solve the optimization problem 
eq. ( p2| ) in each step of the algorithm. The most obvious method at hand to overcome this difficulty is to use Dykstra's 
algorithm [T71[2T]. Dykstra's algorithm is an elaborate technique to compute the projection of X onto the intersection 
of several convex sets by means of modified iterative projections onto the individual sets. To apply this procedure to 
the computation onto C we need to determine the projections onto the set of positive semi-definite Hermitian matrices 
and onto the affine set given by the constraints imposed by the canonical (anti)-commutation relations and other 
symmetries. The former can be implemented straightforwardly. For the latter we were able to determine formulae for 



the submatrices in eq. ( 13 1 and thus gain a method to compute Pc '■ % — > C efficiently. Note that since C is convex, 
Pc (X) is unique for all X G % [T7]. Further, the proposed algorithm finds a lower bound on the ground state energy 
by only considering boundary points of the feasible set and hence is slightly related to the algorithm proposed in [2] . 



There are two parameters defining the accuracy of the estimated lower bound. In principle the algorithm (23) 
converges for perfect projections Pc (X) to a matrix leC minimizing the function E (X) in the limit k —> oo, where 
k is the iteration number. As for all iterative algorithms we need to introduce stopping criteria, on the one hand, for 
the projected gradient algorithm and, on the other hand, for Dykstra's algorithm. For the latter we use the robust 
stopping criterion introduced in [21 j and determined by a small number To y kstra G K.. For the former one can exploit 



the fact that the dual problem of the optimization problem eq. ( 20 ) , which we identify as the primal problem, can be 



solved using the projected gradient algorithm, too. Note that the primal problem can be written as 

minimize cfix 

subject to F + x t F t > 0, 1 j 

where we identify the matrices G and X with the vectors g and x such that E (X) — c = (G, X) = (fx. The Hermitian 
matrices F Q , {F i } i comprise the linear constraints dictated by the fermionic anti-commutator relations such that the 



matrix Fq + ^ XiFi is an element of the affine set for all Xi G K. The dual problem of eq. (24) is then for the dual 



variable Z G U [TO] 

maximize — (Fq,Z) 

subject to gi = (Fi,Z) (25) 
Z > 0. 

Note that the constraints of the dual problem force the solution Z* to lie both in the positive semi-definite cone and 
in the affine set given by the linear equations. The intersection of these two convex sets is again convex and the dual 
problem becomes 

minimize (F ,Z) , , 

subject to Z G £>, 1 > 

where T> denotes the intersection of the positive semi-definite cone with the set given by the linear constraints 
gi = (Fi, Z). Finally, lemma [I] shows that the projected gradient algorithm can be applied to solve the dual problem 
as well. 

The benefit of the formulation as a primal-dual problem is that the solution of the dual problem is a lower bound on 
the solution of the primal problem. Hence, solving the primal and dual problems simultaneously using the projected 
gradient algorithm provides us with an ingenious certificate for the accuracy of the solution. Consequently, one way 
to introduce a stopping criterion for the projected gradient algorithm is to initiate a small quantity Tp r i ma i—d ua i and 
stop the computations when 

E (Xk) — H (Zk) < Tprimal-dual, (27) 



where H (Z) = —(F ,Z) + c is basically the dual function up to a constant. Note that eq. (271 is equivalent to 
the duality gap of the primal and dual problem and that for strictly feasible primal problems, i.e. for optimization 
problems of the form eq. (|24[) for which there exists an x such that Fq + XiFi > 0, the duality gap is zero if 



evaluated at optimal points [10] . For the optimization problem eq. (20), one can easily find a point in the feasible set 



which is strictly feasible. Consequently, we have E (X*) — H (Z*) = 0, where X* and Z* are solutions of the primal 



and dual problems, respectively, and hence eq. (27) suits perfectly as a stopping criterion. Due to the numerical 
cost to solve the primal and dual problems simultaneously, we decide on using a different stopping criterion for the 
projected gradient algorithm. We stop when the improvement per iteration step, i.e. 

E(X k+1 )-E(X k ) <t e , (28) 

falls below a predetermined small number te G K. We implemented both routines for systems considering only second 
moments and numerical experiments suggest that for a sufficiently small number te the latter stopping criterion yields 
results that are as good as in the primal-dual formalism. 

Another degree of freedom in the formulation of the algorithm is the step length a. Numerical examples propose 
that for the first step in the direction of the negative gradient a should be chosen to be large. For the remaining steps 
the computational cost could be reduced by choosing a small but appropriate value of a [22j . 



IV. APPLICATIONS 



In the remainder of this work, we provide a number of numerical examples that demonstrate that the proposed 
method can be applied to rather large systems and is capable of extracting accurate information concerning the 
ground state energy as well as the correlation functions. The purpose of what follows is not to improve lower bounds 
on unknown ground state energies, but rather to analyse the behaviour of the proposed algorithm for finite but large 
lattice models. 



A. Second Moments and the Ising Model 



We begin our exposition of the numerical application of the algorithm with a simple, exactly solvable model for which 
the relaxation of the energy minimization problem has been analysed in great detail exploiting different numerical 
methods in [26] . The Hamiltonian we are considering as a first application of the projected gradient algorithm is given 
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by 



AT-l 
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(29) 



k = l 



where the operators describe spinless fermions. The ground state of this model coincides with that of the Ising model 
in a transverse magnetic field with periodic boundary conditions 
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(30) 



when j c < and AT is even. This can be seen from the connection of eq. (30) with eq. (29 1 via a Jordan- Wigner 
transformation [231 124|. To analyse a translational invariant system with periodic boundary conditions, we simply 
change the signs of the second and the fourth terms in the Hamiltonian eq. (29 1. 



Now we require only constraints that are due to second moments, i.e. we require that 



(31) 



and that the entries fulfil the constrains dictated by the anti-commutator relations. Solving the relaxed minimization 
problem eq. (20) by applying the projected gradient algorithm and comparison of the so-obtained lower bound to the 



exact ground state energy [23l [24] 



yields the results presented in table 
eq. ( 29 ) is computed by the formula 
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that the ground state energy for the system given by the Hamiltonian 
with parameter If. = 2k + 1, while one chooses If. = 2k for the energy of a 



translational invariant system with periodic boundary conditions. 
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TABLE I. Lower bounds on the ground state energies for the translational invariant Jordan- Wigner transformation of the Ising 
Hamiltonian with parameters j c = —1 and h = 1/2 for modes. The Hamiltonian is identical to eq. (29 1 but the signs of 
the second and fourth terms are changed to make the Hamiltonian translational invariant. The relative deviation AE is of the 
order of the parameters of the stopping criteria TDykstra =te = 10~ 8 , which limits the accuracy of this method. We initialized 
the algorithm with a = N for the first step and a <C iV for the remaining steps until convergence. 

In fact, it can be shown directly that the ground state energy of the Ising model is exactly reproduced with the 
present method as the Ising model can be diagonalized by an orthogonal mode transformation mapping second 
moments to second moments, this making the transformed optimization problem straightforward to solve [25] . Thus 
the so obtained lower bounds are, in this case, superior to other approaches such as Anderson lower bounds. In 
particular, solving the optimization problem for the non-translational invariant Hamiltonian eq. (29) for even and 

d30b. Even for 



negative j c yields the exact ground state energies for the translational invariant Ising Hamiltonian eq. 
this Hamiltonian one is able to solve the optimization problem for large particle numbers since one is only considering 
the 1-positivity conditions which scale linearly in the system size. 



B. Fourth Moments and the Heisenberg model 



While the Ising model is quasi-free it is of considerable interest to study interacting models, such as the Heisenberg 
Hamiltonian, that contain, for example, fourth-order fermionic operators (see [27] for the study of another interacting 
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fermionic system where a different numerical approach is employed). In the following, we consider the performance 
of the proposed algorithm in this situation. The Heisenberg Hamiltonian 



N-l 



(33) 



can be transformed to spinless fermions by the Jordan- Wigner transformation. In this picture, the Hamiltonian reads 

N-l N-l 



h = 2 j ^2 [ a l+i a i + a l a i+i [ a !+i ai +! + a I a « 

1=1 t=l 

N-l 

+4J 4 a i+l a i+l a i + J ■ (N -1) . 



(34) 



Now we are in a position to compare the exact ground state energies of the Heisenberg model with the lower bounds 
obtained by applying the suggested algorithm. Indeed, for the Heisenberg Hamiltonian with parameter J = 1/2 the 
lower bounds on the ground state energy lie very close to the actual energies. Table |ll| presents some numerical results. 



N 


E exa ct 


Lower bound 


Ratio 


6 


-4.9872 


-4.9974 


0.9979 


8 


-6.7499 


-6.7739 


0.9964 


10 


-8.5161 


-8.5557 


0.9954 


12 


-10.2842 


-10.3412 


0.9945 


14 


-12.0534 


-12.1291 


0.9938 


16 


-13.8235 


-13.9189 


0.9931 


18 


-15.5940 


-15.7099 


0.9926 


20 


-17.3649 


-17.5021 


0.9922 



TABLE II. Lower bounds for the ground state energy of the Heisenberg Hamiltonian eq. ( |34[ ) of N modes and with parameter 
J = 1/2. The exact ground state energy is obtained by diagonalizing the Hamiltonian in the spin picture. The parameters 
defining the stopping criteria for Dykstra's algorithm and for the projected gradient algorithm are TDykstra = 10 -4 and 



TE 
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Increasing the precision of the algorithm by decreasing these parameters is accompanied by an improvement of the 



lower bounds and an increase of the computing time. 



The range of applicability of the projected gradient algorithm can be extended considerably by analysing a system 
that is translationally invariant. To this end, we make the Hamiltonian in eq. (34) translationally invariant with 
periodic boundary conditions. The results for the ground state energy of this Hamiltonian are presented in the left- 
hand side of figure |] As a reference, we depict the exact ground state energies for this model for up to N = 22. 
Additionally, the right-hand side of figure [2] presents the nearest-neighbour and next-nearest-neighbour correlations 
determined by the gradient projection algorithm. 

To illustrate the rate of convergence of the algorithm for fourth moments, we present the improvement per iteration of 
the function E {X k ) = (G, Xk) + c for the translational invariant Jordan- Wigner transformed Heisenberg Hamiltonian 
describing N = 50 sites in figure [3j It is noticeable that the algorithm converges to a good estimation in only a few 
steps and that the most time-consuming factor is given by the projection onto the set C. 

Note that in all calculations we initialized the algorithm with the zero matrix since we wanted to test the algorithm 
without imposing prior knowledge. Generally, for translational invariant systems it is advantageous to complete the 
circulant submatrices obtained from a run with site number iVi to circulant matrices describing N2 sites and use these 
matrices as an initial matrix for N2 > N±. 

Furthermore, since the matrices comprising the fourth moments in the translational invariant case decompose into 
blocks of diagonal matrices, it is possible to parallelize the projection onto the positive semi-definite matrices. In fact, 
this step in the computation is the most expensive since one has to diagonalize the matrices in each step. Parallelizing 
this computation might speed up the calculations. So far, no effort has been made to do so since our intention is 
rather to demonstrate that this technique yields good results. 
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11.32 
0.30 
0.28 
0.26 
11.24 
0.22 

-0.58 
-0.60 
-0.62 
-0.64 



estimated correlations 
exact correlations 

W<+J>oo 



10 14 IS 22 26 30 34 38 42 46 50 

N 



FIG. 2. The system under consideration is the Hamiltonian eq. (34 1 extended to a translational invariant system with periodic 
boundary conditions and parameter J = 1/2. The stopping criteria are initialized with TDykstra = 10 -4 and te = 10 -4 . 
Left: lower bounds on the ground state energies ejv = En/N per site. Circles indicate the estimated lower bounds obtained 
via the projected gradient algorithm. The exact energy per particle is illustrated as triangles for up to N — 22 particles. 
Additionally, the energy per particle in the thermodynamical limit = — 4|J|(ln2 — 0.25) 23 is shown as a horizontal 
line to indicate where the exact energies converge to. Right: nearest- and next-nearest-neighbour correlations. Shown are 
the estimated correlations as circles, the exact correlations as triangles and the correlations in the thermodynamical limit as 
lines, i.e. (of^i)^ = (1 - 41n2)/3 = -0.59086292 and {tr|(r| ! +2 ) 00 = (1 - 16 In 2 + 9C(3))/3 = 0.24271908 where C(3) = 
| X]fcli( — l) fc_1 /[fc 3 ( 2 ^)\ — 1.2020569.... The nearest-neighbour correlations are obtained by averaging over (a™cr" +1 ) for 



a — x,y,z, which is possible since the system 
proceeded the same way. 



rotationally invariant. 



For the next-nearest-neighbour correlations, we 



3.78 
3.76 



3.72 

3.70 - 

3.68 - 

3.66 Li 



(I 



FIG. 3. Energies for each iteration step for the Hamiltonian eq. (341 extended to a translational invariant system for N = 50 
particles illustrating the good behaviour of the gradient projection algorithm with respect to the iteration step k. Note that 
the data is plotted double logarithmically. 



V. CONCLUSIONS 



In this work, we have presented a method for the determination of lower bounds on the ground state energy of 
quantum many-body systems. To this end, we followed a large body of work to relax the general ground state energy 
finding problem to obtain a semi-definite programme that involves a number of variables that scale polynomially 
with the system size. To overcome unfavourable scaling in time and memory of standard primal-dual interior point 
solvers, we combined a projected gradient algorithm with Dykstra's algorithm to obtain access to larger system sizes 
and proved its convergence. For quasi-free systems the method yields the exact ground state energies, while for 
general interacting models it provides very good lower bounds and direct access to correlation functions between sites. 
Different properties of the underlying system such as higher spatial dimensions jTHj , bosonic degrees of freedom 28i , 
symmetries (29) and translational invariance can be included straightforwardly and the relative weight of higher order 
correlation functions can be adjusted to optimize for computational performance. 

We hope that further investigations will help to optimize this approach and thus elucidate its power and limitations. 
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Appendix A: Proof of Convergence 



The analysis of lemma [l] requires some basic properties of the projection operators defined in eq. (22 1. The following 
results can be found, for example, in [IT], theorems 4.1 and 5.5. 

Lemma 2 Let Pq :% — > C be the ■projection onto the convex set C. 

1. For all Z eC and all X E H we have (P c (X) - X,Z - P c (X)) > 0. 

2. For all X, Y G H we have (P c (X) - P c (Y) ,X — Y) > 0. 

The following proof is motivated by 20} . 
Proof of lemma [T] 

Before we start with the main proof, we need to establish the relation 



1/ £ ~ 

for a > and where the sequence {X k } k is defined via 



<G, X k - X k+1 ) > IW1 ~ k " F (Al) 



X k+1 = P c (X k - aG) , (A2) 
i.e. the projected gradient algorithm. To verify relation ( |Al[ ), note that with the first inequality of lemma [2j we find 

(^fe+i - Xk + ctG, X k - X k +i) 

= (P c (X k - aG) -X k + aG, X k - P c (X k - aG)) (A3) 



> 



since X k £ C. Hence we have 



(aG,X k ~ X k+1 ) >\\X k - X k+1 \\% (A4) 
and relation (Al) is established. Note that with this expression the following holds for all a > 0: 

(G,X k+1 -X k ) <0. (A5) 
In a first step, we show that the sequence E (X k ) is monotonically decreasing. We find that 

(G,X k+1 ) < (G,X k ), 

^(G,X k+1 )+c<(G 1 X k )+c, (A6) 
^E(X k+1 ) <E(X k ), 

where E (X) = (G, X) + c, which holds for all a > 0. Hence the sequence E (X k ) is monotonically decreasing. Since 
we assume that the function E (X) is lower bounded on C, we then conclude that this sequence converges for all 
a > 0, i.e. E (X k ) — > E* for k — > oo. It remains to be shown that E* is indeed the minimum of the function in C. 
In a second step, we show that 

lim||X fe+1 -X fe || F = 0. (A7) 

k— >oo 



Since E (X) is linear in X, we have 



E (X k ) - (G, X k )+c+ (G, X k+1 ) - (G, X k+1 ) 

= (G,X k+l )+c+(G,X k -X k+l ) (A8) 
= E (X k+ i) + (G, X k — X k+ i) 
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and hence 



limE(X k )= UmE(X k+1 ) + lim (G,X k - X k+1 ). (A9) 

k— »oo k— >oo k^-oo 



Since the sequence E (X k ) converges, i.e. E (X k ) — > E* , we find that 



lim(G,X k -X k+1 )=0. (A10) 

k—±oo 



Inequality (Al) then yields 



lim (G,X k - X k+1 ) > lim l|Xfc+1 ^ ., 

k— >oo k— Yoo (X 

^ > lim l|Xfe+1 ~ Xfc|1 ^ 

fe— >oo a 



(All) 



and hence 



Thus we have 



\\m\\X k+1 -X k \\ 2 F = Q. (A12) 

k— voo 



lim||A fe+1 -X k \\ F = 0. (A13) 

In a third step, we show that the elements X k are bounded. Let X £ C be such that E (Xj < E (X) for all leC. 
Then 

ll-^As ~ ^Hf = ll-^fc — + Afe + i — 

= \\X k — ^fe+l||i? + H-Xfc+l — -X"! [ jr 
+2(Afc — Xfe+i, Afc + i — A) 

= ||A/c — Xfc+iUjr + | — X\\ F 

+2(X k -aG + aG- X k+1 ,X k+1 - X) 
= \\X k -X k+1 \\ 2 F + \\X k+l -Xf F (A14) 
+2(X k -aG-Pc (X k - aG) , P c {X k - aG) - X) 
+2(aG 1 P c (X k -aG)-X) 

= ||Afe — -Xfc+iUjr + ||^fe+l — X\\ F 

+2(P C (X k - aG) - X k + aG, X - P c (X k - aG)) 
+2a [E (X k+1 )-E(X)] , 

since (G, P c (X k - aG) - X) = (G, X k+X - X) = E (X k+ i) - E (A) . Note that with lemma [2] and the fact that 
E (A) < E (A) for all X 6 C, one can show that the last two terms are greater than or equal to zero. Therefore 

\\X k -X\\ F >\\X k -X k+1 \\ 2 F + \\X k+1 -X\\ F (A15) 

and consequently 

||A fc+1 - X\\ 2 F < \\X k - X\\% -\\X k - X k+1 \\ 2 F , (A16) 

which leads to 

\\X k+1 -X\\ F <\\X k -X\\ F . (A17) 

We conclude that the elements X k are bounded. As a consequence we know that, since C is closed, there is at least 
one converging subsequence A„ fe — ¥ X* for k — ¥ 00, where A* is an accumulation point of the sequence {X k } k . 

In a last step, we show that for any accumulation point X* it holds that ^(A*) < E (X) for all X E C. With 
lemma [2] we find for all Z € C 

(X k+1 -X k + aG,Z-X k+1 )>0, 

•4=> a(G, X k+ i — Z) < (X k+ i - X k , Z - X k+ i). 
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Further, we have 



oc(G, Xk+i 


-z) 








< (Xi , i - 


Xi z - 

^kj 


Xi, 

^ l k+ 


1/ 




< (x k+1 - 


Xk,Z — 


x k - 


- Xk — Xk- 


hi) 


< (x k+1 - 


Xk,Z — 


x k ) 


+ (Xk+i - 


Xk,X 


< (Xk+i — 


Xk,Z — 


x k ) 


- \ \Xk+i - 


Xk\\p 


< (Xk+i — 


Xk,Z — 


x k ) 






< Afc + 1 - 


- Xk\\F ■ 


\\z - 


Xk \\f> 





(A19) 



where we exploited the Cauchy-Schwarz inequality. Finally, we find that 



(G, X k — Z) — (G, X k — X k +i) + (G, X k +i — Z) 

< (G,X k - Xk+i) + -\\X k +i - X k \\ F -\\Z- X k \\ F 
a 



(A20) 



holds for all Z G C. Now choose a converging subsequence X nk such that X„ 



X* . We already know that 



■ X„ 



for k — > oo. 



(G,X nk — X nk +i) —> due to equation (A10). Further, equation (A13) claims that HA^+i 
Additionally, since the sequence {Xk} k is bounded, we know that for X £ C, where E (X) < E (X) for all X £ C, 



x nk \\ F = \\z-x + x-x nh \\ F 

<\\Z-X\\ F + \\X-X nk \\ F 

< \\Z - X\\ F + \\X - X \\f 
<C(Z), 



(A21) 



where we exploited that the sequence | \X — X nk \\ F is upper bounded by \\X — A | \ F , see equation ( A17 1 , and where 
C (Z) is a constant which may depend on arbitrary Z £ C. Hence we find for all Z £ C and all accumulation points 

X* 

lim (G,X nk -Z) <0, 

k— >oo 

=> (G, X* -Z)<0, 

=> (G, X*) < (G, Z), (A22) 
=► (G, X*) + c < (G,Z)+c, 
E(X*) <E{Z). 

Since we know that the sequence E (Xk) converges, we can choose the subsequence defined by X nk — > X* to see that 
E (X nk ) — > E* = E(X*). With the argumentation above, we know that E* < E (Z) for all Z £ C and hence the 
algorithm ( A2) finally ends up at a point minimizing the function E (X) in the convex set C. □ 



Appendix B: Constraints for a Fermionic System for up to 2-positivity Conditions 



Here we list the constraints on the matrices for the second and fourth moments which arise by applying the 
fermionic anti-commutator relations and further specify the matrices under the additional assumption of particle 
number conservation for a system composed of two spinless fermions. For the second moments of a general, fermionic 
system, we find 

Tkd = (a k ai), 

Sk,i = (a k a\)=5 k>l -T hk , (Bl) 
U k ,i = (a k ai) = -Uim- 
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The fourth moments yield 



Gkl,mn 
Hkl.mn 



{ a k a j a n a m) — —Mik,mn — —^kl,t 



— kn.inl — 


krfi.ln 


— k n,lm 5 




= {akaia n am} 










-^-km,nl ~ 


*i km, In 


^kn,lm 




-H-nl,mk 




H rim, 


— ^nl,km ~~ 


Hyik^lm — 


H n7T i,lk — 


-H-lk ,rnn 


— tllk.nrn ~~ 


-H-lmjik — 


771, kn — 


-Eiln,km 


— ^ln,mk ~~ 


Hnil,kri — 


Hml,nk ~ 


ffmk,7il 


-H-mk,ln ~~ 


"mTi,^ — 


"mn,fcl ) 





kl 



(B2) 



.mn 
*kl,mn 
Qkl .mn 



(alaial l a m ) = 5^ n T k ^ r 
(akcualdm) = 6i <n U kt , 



Ui m + G n k,mi and 



(akaialaln) 



+Si m T n .k — fik, m T n .i + M n 



Ik- 



Assuming particle number conservation for a system composed of two spinless fermions, the matrices containing 
the second moments read 



T 



( ( a i a i) ( a l a 2)\ 
\(a\a 2 )* {a\a 2 )J 



\-{a\ ai ) -{a\a 2 )*\ 



-(a[a 2 ) 1 



{<Aa 2 )J ' 



while the matrices for the fourth moments are 



M 







a\a\a 2 ai) — (a\a\a 2 ai) 



t„t. 



— {a\a\a 2 a\) (a\a\a 2 a-\) 



t„t. 







0/ 



/ (a[ai) 


(a|a 2 )* 
\(a\a\a 2 ai) 



(a\ai) - (a{a 2 a2ai) 


(a\a 2 )* 



(a\a 2 ) 


(0^02) - (a|a 2 a 2 ai) 




(a\ala 2 ai)\ 
(a\a 2 ) 


(a 2 a 2 ) / 



(B3) 



/o 0\ 

1 + (a|a 2 a 2 ai) - (ajai) - (a 2 a 2 ) -1 - (a|a 2 a 2 ai) + (ajai) + (a 2 a 2 ) 

—1 — (aja^c^ai) + (a{ai) + (a 2 a 2 ) 1 + (aja^^a!) — (a\ai) — (a\a 2 ) 

\o " " 0/ 

The presented matrices T, S, M, R and Q now satisfy the fermionic anti-commutator relations for two spinless fermions. 



For a system where the particle number is conserved the constraint on the relaxed optimization problem eq. (12 1 is 



that these matrices are all simultaneously positive semi-definite, which generates the 1- and 2-positivity conditions. 
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